clc;
clear;

%% 设置参数
dim = 2; % 维度(1 or 2)
interp_K = 4; % 插值次数
interp_type = "P"; % 插值类型(P or Q)仅二维有效




%% 计算部分
if dim == 1
    left = -1;
    right = 1;
elseif dim == 2 && interp_type == "Q"
    left = -1;
    right = 1;
elseif dim == 2 && interp_type == "P"
    left = 0;
    right = 1;
end
K = interp_K;
reference_Point_1D = linspace(left, right, K+1);
if dim == 1
    % (-1,1)区间K次插值坐标
    fprintf("参考单元[%d,%d]上%d次插值坐标:\n", left, right, K);
    disp(reference_Point_1D);
    fprintf("插值点见图\n");
    figure("WindowStyle", "docked");
    scatter(reference_Point_1D, zeros(size(reference_Point_1D,2)), "r*");
    text(reference_Point_1D', zeros(size(reference_Point_1D,2),1), int2str((1:size(reference_Point_1D,2))'), "FontSize", 10, "FontWeight", "bold", "Color", "k");
    axis("equal", "tight");
    % 一维K次多项式
    str1 = '1';
    for i = 1:K
        str1 = [str1 ' + x^' num2str(i)];
    end
    fprintf("一维%d次多项式:\n", K);
    disp(str1);
    % 各阶导数
    % syms x
    % P = sym(str2sym(str1));
    % for i = 1:K
    %     fprintf("其%d阶导数:\n",i);
    %     disp(diff(P,x,i));
    % end
    % 基函数矩阵(求K个基函数)
    str_split = strsplit(strrep(str1, ' +', ''), ' ');
    str_A = cell(K+1, K+1);
    for i = 1:K+1
        for j = 1:K+1
            str_A{i, j} = str_split{1, j};
        end
    end
    syms x
    A = sym(str2sym(str_A));
    for i = 1:K+1
        for j = 1:K+1
            A_(i, j) = subs(A(i, j), 'x', reference_Point_1D(i));
        end
    end
    % fprintf("%d次插值基函数的矩阵:\n", K);
    % disp(A_);
    fprintf("其逆矩阵(即各项系数):\n");
    A_inv = inv(A_);
    disp(A_inv);
    t1 = sym2cell(A_inv);
    t2 = cell(K+1,K+1);
    for j = 1:K+1
        for i = 1:K+1
            t2{i,j} = t1{i,j}*str2sym(str_split(1,i));
        end
    end
    sym_Phi = cell(K+1,1);
    for i = 1:K+1
        u = sym(0);
        for j = 1:K+1
            u = u + t2{j,i};
        end
        sym_Phi{i,1} = u;
    end
    fprintf("一维%d次基函数:\n", K);
    disp(sym_Phi);
elseif dim == 2
    if interp_type == "P"
        ND_tri = 0.5*(K+1)*(K+2);
        % (0,1)x(0,1)三角形区域k次插值坐标
        reference_Point_2D_Tri = zeros(ND_tri,2);
        j = 1;
        k = 1;
        for i = K+1:-1:1
            reference_Point_2D_Tri(j:j+i-1,1) = reference_Point_1D(1:i);
            reference_Point_2D_Tri(j:j+i-1,2) = reference_Point_1D(k);
            j = j + i;
            k = k + 1;
        end
        fprintf("三角形参考单元P_%d次插值坐标:\n", K);
        disp(reference_Point_2D_Tri);
        fprintf("插值点见图\n");
        figure("WindowStyle", "docked");
        scatter(reference_Point_2D_Tri(:,1), reference_Point_2D_Tri(:,2), "r*");
        text(reference_Point_2D_Tri(:,1), reference_Point_2D_Tri(:,2), int2str((1:size(reference_Point_2D_Tri,1))'), "FontSize", 10, "FontWeight", "bold", "Color", "k");
        axis("equal", "tight");
        % 二维P_k次多项式
        str2 = '1';
        for i = 1:K
            str2 = [str2 ' + x^' num2str(i)];
        end
        for j = 1:K
            str2 = [str2 ' + y^' num2str(j)];
        end
        for i = 1:K-1
            for j = 1:K-1
                if i + j <= K
                    str2 = [str2 ' + x^' num2str(i) '*y^' num2str(j)];
                end
            end
        end
        fprintf("二维P_%d次多项式:\n", K);
        disp(str2);
        % 各阶导数
        % syms x y
        % Pk = sym(str2sym(str2));
        % for j = 1:K
        %     fprintf("其偏x的%d阶导数:\n",j);
        %     disp(diff(Pk, x, j));
        %     fprintf("其偏y的%d阶导数:\n",j);
        %     disp(diff(Pk, y, j));
        % end
        % 基函数矩阵(求K个基函数)
        str_2_split = strsplit(strrep(str2, '+ ', ''), ' ');
        str_Tri = cell(ND_tri, ND_tri);
        for i = 1:ND_tri
            for j = 1:ND_tri
                str_Tri{i, j} = str_2_split{1, j};
            end
        end
        syms x y
        A_tri = sym(str2sym(str_Tri));
        k = 1;
        for i = 1:ND_tri
            for j = 1:ND_tri
                A_tri_(i,j) = subs(A_tri(i,j), 'x', reference_Point_2D_Tri(k,1));
                A_tri_(i,j) = subs(A_tri_(i,j), 'y', reference_Point_2D_Tri(k,2));
            end
            k = k + 1;
        end
        % fprintf("三角形区域参考单元%d次插值基函数的矩阵:\n", K);
        % disp(A_tri_);
        fprintf("其逆矩阵(即各项系数):\n");
        A_inv_tri = inv(A_tri_);
        disp(A_inv_tri);
        t1 = sym2cell(A_inv_tri);
        t2 = cell(ND_tri, ND_tri);
        for j = 1:ND_tri
            for i = 1:ND_tri
                t2{i,j} = t1{i,j}*str2sym(str_2_split(1,i));
            end
        end
        sym_Phi_tri = cell(ND_tri,1);
        for i = 1:ND_tri
            u = sym(0);
            for j = 1:ND_tri
                u = u + t2{j,i};
            end
            sym_Phi_tri{i,1} = u;
        end
        fprintf("二维P_%d次Phi函数:\n", K);
        disp(sym_Phi_tri);
    elseif interp_type == "Q"
        ND_rec = (K+1)*(K+1);
        % (-1,1)x(-1,1)矩形区域k次插值坐标
        reference_Point_2D_Rec = zeros(ND_rec,2);
        j = 1;
        for i = 1:K+1:ND_rec
            reference_Point_2D_Rec(i:i+K,1) = reference_Point_1D(:);
            reference_Point_2D_Rec(i:i+K,2) = reference_Point_1D(j);
            j = j + 1;
        end
        fprintf("矩形单元参考单元Q_%d次插值坐标:\n", K);
        disp(reference_Point_2D_Rec);
        fprintf("插值点见图\n");
        figure("WindowStyle", "docked");
        scatter(reference_Point_2D_Rec(:,1), reference_Point_2D_Rec(:,2), "r*");
        text(reference_Point_2D_Rec(:,1), reference_Point_2D_Rec(:,2), int2str((1:size(reference_Point_2D_Rec,1))'), "FontSize", 10, "FontWeight", "bold", "Color", "k");
        axis("equal", "tight");
        % 二维Q_k次多项式
        str3 = '1';
        for i = 1:K
            str3 = [str3 ' + x^' num2str(i)];
        end
        for j = 1:K
            str3 = [str3 ' + y^' num2str(j)];
        end
        for i = 1:K
            for j = 1:K
                str3 = [str3 ' + x^' num2str(i) '*y^' num2str(j)];
            end
        end
        fprintf("二维Q_%d次多项式:\n", K);
        disp(str3);
        % 各阶导数
        % syms x y
        % Qk = sym(str2sym(str3));
        % for j = 1:K
        %     fprintf("其偏x的%d阶导数:\n",j);
        %     disp(diff(Qk, x, j));
        %     fprintf("其偏y的%d阶导数:\n",j);
        %     disp(diff(Qk, y, j));
        % end
        % 基函数矩阵(求K个基函数)
        str_3_split = strsplit(strrep(str3, '+ ', ''), ' ');
        str_Rec = cell(ND_rec, ND_rec);
        for i = 1:ND_rec
            for j = 1:ND_rec
                str_Rec{i,j} = str_3_split{1,j};
            end
        end
        syms x y
        A_rec = sym(str2sym(str_Rec));
        k = 1;
        for i = 1:ND_rec
            for j = 1:ND_rec
                A_rec_(i,j) = subs(A_rec(i,j), 'x', reference_Point_2D_Rec(k,1));
                A_rec_(i,j) = subs(A_rec_(i,j), 'y', reference_Point_2D_Rec(k,2));
            end
            k = k + 1;
        end
        % fprintf("矩形单元参考单元%d次插值基函数的矩阵:\n", K);
        % disp(A_rec_);
        fprintf("其逆矩阵(即各项系数):\n");
        A_inv_rec = inv(A_rec_);
        disp(A_inv_rec);
        t1 = sym2cell(A_inv_rec);
        t2 = cell(ND_rec, ND_rec);
        for j = 1:ND_rec
            for i = 1:ND_rec
                t2{i,j} = t1{i,j}*str2sym(str_3_split(1,i));
            end
        end
        sym_Phi_rec = cell(ND_rec,1);
        for i = 1:ND_rec
            u = sym(0);
            for j = 1:ND_rec
                u = u + t2{j,i};
            end
            sym_Phi_rec{i,1} = u;
        end
        fprintf("二维Q_%d次Phi函数:\n", K);
        disp(sym_Phi_rec);
    end
end